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SYNCHRONISING AND NON-SYNCHRONISING DYNAMICS FOR A 
TWO-SPECIES AGGREGATION MODEL 

CASIMIR EMAKO-KAZIANOU, JIE LIAO, AND NICOLAS VAUCHELET 


Abstract. This paper deals with analysis and numerical simulations of a one-dimensional two- 
species hyperbolic aggregation model. This model is formed by a system of transport equations 
with nonlocal velocities, which describes the aggregate dynamics of a two-species population in 
interaction appearing for instance in bacterial chemotaxis. Blow-up of classical solutions occurs 
in finite time. This raises the question to define measure-valued solutions for this system. To this 
aim, we use the duality method developed for transport equations with discontinuous velocity to 
prove the existence and uniqueness of measure-valued solutions. The proof relies on a stability 
result. In addition, this approach allows to study the hyperbolic limit of a kinetic chemotaxis 
model. Moreover, we propose a finite volume numerical scheme whose convergence towards 
measure-valued solutions is proved. It allows for numerical simulations capturing the behaviour 
after blow up. Finally, numerical simulations illustrate the complex dynamics of aggregates until 
the formation of a single aggregate: after blow-up of classical solutions, aggregates of different 
species are synchronising or nonsynchronising when collide, that is move together or separately, 
depending on the parameters of the model and masses of species involved. 


1. Introduction 

Aggregation phenomena for a population of individuals interacting through an interaction 
potential are usually modelled by the so-called aggregation equation which is a nonlocal nonlinear 
conservation equation. This equation governs the dynamics of the density of individuals subject 
to an interaction potential K. In this work, we are interested in the case where the population 
consists of two species which respond to the interaction potential in different ways. In the 
one-dimensional case, the system of equations writes: 

(1.1) dtPa + Xadx {a{p)pa) = 0, for 0 = 1,2, 

with 

a{p):= / dxK{x - y)p{t,dy), p := Oipi + e 2 P 2 , 

Jr 

where Oa,Xa for a = 1,2 are positive constants. 

In this work, we are interested in the case where the interaction potential K in (jl.ip is pointy 
i.e. satisfies the following assumptions: 

(HI) KeC\M\{0}). 

(H2) Vx G M, K{x) = K{-x). 

(H3) dxK^L^{R). 

(H4) K is A-concave with A > 0 i.e., 

Vx,y G M*, {dxK{x) - dxK{y)) {x - y) < X{x - yf. 

The aggregation equation arises in several applications in biology and physics. In fact, it is 
encountered in the modelling of cells which move in response to chemical cues. The velocity of 
cells a{p) depending on the distribution of nearby cells represents the gradient of the chemical 
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substance which triggers the motion. Cells gather and form accumulations near regions more 
exposed to oxygen as observed in |20l [M] . We can also describe the movement of pedestrians 
using the aggregation equation as in [16] where the velocity of pedestrians is influenced by the 
distribution of neighbours. This equation can also be applied to model opinion formation (see 
|25] 1 where interactions between different opinions can be expressed by a convolution with the 
kernel K. 

From the mathematical point of view, it is known that solutions to the aggregation equation 
with a pointy potential blow up in finite time (see e.g mum)- Then global-in-time existence 
for weak measure solutions has been investigated. In [5|, existence of weak solutions for single 
species model has been obtained as a gradient flow. This technique has been extended to the 
two-species model at hand in m- Another approach of defining weak solution for such kind of 
model has been proposed in mm for the single species case. In this approach, the aggregation 
equation is seen as a transport equation with a discontinuous velocity a{p). Then solutions in 
the sense of duality have been defined for the aggregation equation. 

Duality solutions has been introduced for linear transport equations with discontinuous ve¬ 
locity in the one-dimensional space in [3]. Then it has been adapted to the study of nonlinear 
transport equations in Iiii8i[i7|. In [min], the authors use this notion of duality solutions 
for the one-species aggregation equation. Such solutions are constructed by approximating the 
problem with particles, i.e. looking for a solution given by a finite sum of Dirac delta functions. 
Particles attract themselves through the interacting potential K, when two particles collide, 
they stick to form a bigger particle. 

In this work, we extend this approach to the two species case. To do so, we need to modify 
the strategy to the problem at hand. Indeed, collisions between particles of different species 
are more complex: particles can move together or separately after collision. This synchronising 
or non-synchronising dynamics implies several difficulties for the treatment of the dynamics 
of particles. In fact, particles of different species can not stick when they collide. Then an 
approximate problem is constructed by considering the transport equation with the a regularized 
velocity. Then measure valued solutions are constructed by using a stability result. 

An important advantage of this approach is that it allows to prove convergence of finite vol¬ 
ume schemes. Numerical simulations of the aggregation equation for the one-species case, which 
corresponds to the particular case of (HTD when setting p 2 = 0, have been investigated by several 
authors. In [8| the authors propose a finite volume method consistent with the gradient flow 
structure of the equation, but no convergence result has been obtained. In |9], a Lagrangian 
method is proposed (see also the recent work m)- For the dynamics after blow up, a finite vol¬ 
ume scheme which converges to the theoretical solution is proposed in mu- In the two-species 
case, the behaviour is more complex since the interaction between the two species can occur 
and they may synchronise or not i.e. move together or separately depending on the parameters 
of the models and the masses of species. A numerical scheme illustrating this interesting syn¬ 
chronising or non-synchronising dynamics is provided in Section 6. In addition, a theoretical 
result on the convergence of the numerical approximation obtained with our numerical scheme 
towards the duality solution is given. Such complex interactions phenomena have been observed 
experimentally in m- 

System (|I.lh can be derived from a hyperbolic limit of a kinetic chemotaxis model. In the 
case of two-velocities and in one space dimension, the kinetic chemotaxis model is given by 
( 1 . 2 ) 

I dtf^ + vd,,f^ = ^ J^{Ta[S]{v',v)f^{v') -Tc,[S]{v,v')f^{v)) dv', 

\ - + 8 ^ = 01 imi) + fii-i))+ 02 (/i(i) + /K-i)), 


a = 1,2, v€V = {±1}, 
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where f^{x,v,t) stands for the distribution function of a-th species at time t, position x and 
velocity v, S^(t,x) is the concentration of the chemical substance, Ta[S]{v,v') is the tumbling 
kernel from direction v gV to direction v' gV and e > 0 is a small parameter. This tumbling 
kernel being affected by the gradient of the chemoattractant, is chosen as in |12] 

(1.3) Ta[S]{v,v') =ll}a{'^ + XaVdxS) , 

where V’a is a positive constant called the natural tumbling kernel and Xa is the chemosensitivity 
to the chemical S. This kinetic model for chemotaxis has been introduced in m to model the 
run-and-tumble process. Existence of solutions to this two species kinetic system has been 
studied in [Ill- 

Summing and substracting equations (|1.2p with respect to u = ±1 for yields 

(1.4) dtp% + dxJa = 0, 

(1-5) dtJl + dxPa = ‘^ixadxS^pl - J^), a = 1,2, 

where p% := /^(l) + /„(—1) and ;= (/^(l) — /„(—!))• Taking formally the limit e ^ 0 in 
we deduce that ^ XadxS^p% in the sense of distributions. Injecting in (jl.4|) . we deduce 
formally that p^ satisfies the limiting equation: 

( 1 - 6 ) + Xadx{{dxS^)pl) = 0 , 

where satisfies the elliptic equation: 

-dxxS^ + s^ = eip\ + e2pl 

This latter equation can be solved explicitly on M and is given by 
(1.7) 5° = i^*(0iP? + 02P^), K = 

Then we recover system This formal computation can be made rigorous. The rigorous 

derivation of ()1.6p from system (jl.2p will be proved in this work. 

The paper is organized as follows. We first recall some basic notations and notions about the 
duality solutions and state our main results. Section 3 is devoted to the derivation of the macro¬ 
scopic velocity used to define properly the product a{p)pa and duality solutions. Existence and 
uniqueness of duality solutions are proved in Section 4, as well as its equivalence to gradient flow 
solutions. The convergence of the kinetic model (11.21) as e —)■ 0 towards the aggregation model 
dUD-drzi) is shown in Section 5. Finally, a numerical scheme that captures the synchronising 
and non-synchronising behaviour of the aggregate equation is studied in Section 6, as well as 
several numerical examples showing the synchronising and non-synchronising dynamics. 

2. Notations and main results 

2.1. Notations. We will make use of the following notations. Let T > 0, we denote 

• L)|_(M) is the space of nonnegative functions of L^(M). 

• Co(K) is the space of continuous functions that vanish at infinity. 

• -^ioc(I^) is the set of local Borel measures, A4b(M) those whose total variation is finite: 

A4;,(M) = {p G A4ioc(IR), \p\ (K) < Too}. 

• Sm = C([0, T], A4b(K) — cr(A4b(M), Co(M))) is the space of time-continuous bounded 
Borel measures endowed with the weak topology. 
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• 7^2(1^) is the Wasserstein space of order 2: 

7^2(1^) = r 7^ nonnegative borel measures in M s.t |/i| (M) = 1, / |x|^/r(dx) 

I Jm 

• For H G C'(M\{0}), we define H: 

H{x), for X / 0, 


< 00 


H = 


0, else. 


We notice that if K satisfies [(H2)| and |(H4)[ we have by taking y = —x in |(H4)| and using [(H^ 
that, Vx € M, dxK{x)x < Ax^. We deduce that (H4) holds for dxK i.e.: 

(2.1) Vx,?/€M, {dxK{x) - dxK{y)){x-y) <\{x- yf. 

We recall a compactness result on 7\db(K) — cT(7\db(K), C'o(M)). If there exists a sequence of 
bounded measures G 7\4b(M) such that their total variations |/r”| (M) are uniformly bounded, 
then there exists a subsequence of yP' that converges weakly to /i in A4b( 


(2.4) 


2.2. Duality solutions. For the sake of completeness, we recall the notion of duality solu¬ 
tions which has been introduced in [3] for one dimensional linear scalar conservation law with 
discontinuous coefficients. Let us then consider the linear conservation equation: 

(2.2) dtp +dx{b{t,x)p) = in]0,T[xM, 

with T > 0. We assume weak regularity of the velocity field h G L°°(]0, r[xM) and b satisfies 
the so-called one-sided Lipschitz (OSL) condition: 

(2.3) < 7(7)) 7 € L^(]0, T[), in the sense of distributions. 

In order to define duality solutions, we introduce the related backward problem 

J dtp + bdxP = 0, 

\ P{T, ■) =p^ & Lipioci^). 

We define the set of exceptional solutions £ as follows 

£” := {p G Lipic,c{]0,T[x'R) solution to (12.4|) with = O} . 

Definition 2.1 [Reversible solutions to (|2.4jl l. We say that p is a reversible solution to (j2.4p if 
and only if p G Lip;oc(]0,T[xM) satisfies (|2.4D and is locally constant on Ve, where Ve is defined 
by 

Ve := {(7,x) G]0,r[x]R; 3pe G £pe{t,x) 7 ^ 0} . 

Definition 2.2 [Duality solutions to (|2.2I) . see m)- We say p G Sm is a duality solution to 
()2.2p in ]0, T[ if for any 0 < r < T, and any p reversible solution to (|2.4I) compactly supported 

the function p[t,x)p[t,dx) is constant on [0,r]. 

JR 


m X, 


The following result shows existence and weak stability for duality solutions provided that 
the velocity field satisfied the one-sided-Lipschitz condition. 

Theorem 2.3 (Theorem 2.1 in [3]). 

(1) Given p*”® G jVIfe(M). Under the assumption (12.3p . there exists a unique p G Sm, duality 
solution to (12.2p . such that p[0, ■) = p*™. 

(2) There exists a bounded Borel function b, called universal representative of b such that 
b = b a. e., and for any duality solution p, 

dtp + dx[bp) = 0, in the distributional sense. 
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(3) Let {bn)nen be a bounded sequence in L°°(]0, T[xM), with bn ^ b in — w*. 

Assume that dxbn < 7 ”(t), where ( 7 ”)nGN is bounded in L^(]0,T[). Consider a sequence 
Pn G <Sm of duality solutions to 

dtPn + dxibnPn) = 0, in ]0,T[x]R, 

such that Pn{0, •) is bounded in and Pn{0, •) ^ /o*™ € Then pn ^ p in 

Sm, where p G Sm is the duality solution to 

dtp + dxibp) = 0, in ] 0 ,r[xM, p{0,-)=p^"h 

Moreover, bnPn bp weakly in r[xM). 

2.3. Main results. Up to a rescaling, we can assume without loss of generality that the total 
mass of each species is normalized to 1. Then we will work in the space of probabilities for 
densities. 

The first theorem states the existence and uniqueness of duality solutions for system (HI) 
and its equivalence with the gradient flow solution considered in m- 

Definition 2.4. (Duality solutions for system (ll.ip l We say that {pi,p 2 ) G T], 

is a duality solution to (inDif there exists a{p) G L°°((0,T) x R) and 7 G Lj^^{[0,T]) satisfying 

dxh < 7 in the sense of distributions, such that for all 0 < ti < t 2 < T, 

(2.5) dtPa + Xadx{d{p)pa) = 0, for 0 = 1,2, p = Oipi + B 2 P 2 , 

in the sense of duality on (ti,t 2 ) and a{p) = dxK * p a.e. We emphasize that the final datum 
for (12.5p should be t 2 instead of T. 


Then, we have the following existence and uniqueness result: 


Theorem 2.5 (Existence, uniqueness of duality solution and equivalence to gradient flow so¬ 
lution). Let T > 0 and {pff'’, pff'') G 7^2(R)^- Under assumptions (HI)- (Hf), there exists a 
unique duality solution {pi,p 2 ) G (^([O,T], 7 ^ 2 ( 1 ^)^) 

{pi)P 2 ){t = 0) = {pff''■, pff'') such that 


to dni) in the sense of Definition \2.f \ with 
hat 

(2.6) a{p) := [ dxK{x - y)p{t,dy), p = 9ipi + 62 P 2 - 

JR 

This duality solution is equivalent to the gradient flow solution defined in ra¬ 


in our second main result, we prove the convergence of the kinetic model ()1.2p towards the 
aggregation model. 

Theorem 2.6 (Hydrodynamical limit of the kinetic model). Assume that Xafli + 82 ) < 1 for 
a = 1,2. Let T > 0 and {fa,S^) be a solution to the kinetic-elliptic equation (II.2h such that 

/£(t = 0) = ffl^^ and G n L^(R) and [ x‘^f(fl^dx < 00. 

JR 

Then, as e —)■ 0, {ff^,S^) converges in the following sense: 

Pa ■= /a(l) + /«(-!) ^ Pa Weakly in Sm, for a = 1,2, 

^ S in C{[0,T],W^'^{R))-weak, 

where pa is the unique duality solution of (USD and S = K*( 6 ipi-\- 62 P 2 ) given in Theorem \2.5[ 


The condition Xafli + ^ 2 ) < 1 in the previous theorem is needed to guarantee that the 
tumbling kernel TaPS"] defined in (II.3p is positive. 

To conclude this Section on the main results, we emphasize that, a finite volume scheme to 
simulate ()2.5p is proposed in Section 6 and its convergence towards duality solutions is stated 
in Theorem 16.31 
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3. Macroscopic velocity 

In this section, we find the representative a of a for which existence and uniqueness of duality 
solutions hold. To this end, we consider the similar system of transport equations to 
associated to the velocity a” which converges to a. Next, the limit of the product a'^{Pa)p^ for 
a, /3 = 1, 2 is computed. 


3.1. Regularisation. We build a sequence (a"')neN which converges to a by considering the 
sequence of regularised kernels dxK^ approaching dxK. 

Lemma 3.1. Let {dxK^)n£N be the sequence of regular kernels defined by 

1 


dxK'^ix) = 


dxK{x), for |xl > 


n 


ndxK ( - ) X, 
n 


else. 


Then 

and 


drK^ € a 


Mx € M, dxK^{-x) = -dxK^{x), 


< ||<9a;ifi||ooj dxxK^ < A in the distributional sense. 


Proof. From [(Hlj| dxK € C'‘^(M\{0}) and S1I1C0 dxl^ IS COIltlllUOllS 3/tj dz ^, AV0 COIlclu.cl0 
dxK^ G C ^{R). From [(H2) we deduce that dxK is an odd function. Using the definition of 
dxK'^ and |(H3)1 we get that ||cIa;ii'"'||oo < ||i92;ifr||oo- From the construction of dxK"^, we have 
that dxK^ = dxK outside the interval and from |(H4)] one sees dxxK'^ < A in ]R\(—A, A) 

in the sense of distributions. In addition, if we take x = — ^ and y = ^ in|(H4)[ we have that 


ndxK 


n 


< A. 


Since dxxK"' = ndxK (i) in [—^, ^], we conclude that dxxK^ < A in in the sense of 

distributions. Finally, we obtain that dxxK^ < A in the sense of distributions. □ 

In the rest of the paper, the notation dxK^ will refer to the regularised kernels of Lemma [3.II 
Given dxK^, the velocity a"' is defined similarly to (12.6p as 


(3.1) 


ypeSM, a'^ip) ■= [ dxK'^ix -y)p{t,dy). 

Jr 


In the following lemma, we show that if p” and p^ admit weak limits pa and pp respectively 
in Sm^ then the limit of the product a'^{p^)p^ is d{pa)pp. Contrary to [22] where the two- 
dimensional case is considered, this limiting measure does not charge the diagonal. 

Lemma 3.2. For a = 1, 2, let {p^} G Sm be a sequence such that Vn G N, Vt G [0, T], Ip^l {t, M) = 
Mq. Suppose that there exists pa in Sm such that 

Pa Pa weakly in Sm-, 

Then, we have 

^ d{pa)pp weakly in Alb([0, T] x M), a, /3 = 1, 2, 
where a"(-) and a(-) are defined in (l^.dTlD respectively. That is for G Co([0,r] x M), 

[ f dxK'^ix - y)p^{t,dx)pa{t,dy)(j){t,x)dt ^ [ f dxK{x - y)pa{t,dx)pp{t,dy)4){t,x)dt. 
Jo JR2 Jq J^2 
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Proof. Before starting the proof of the lemma, we first introduce some notations which simplify 
the computations 


^2 2 ) P^{t,dx,dy) -.= := | (a^, 2/) G x / y, |x - y| < ^ 

dx, dy) := pa{t, dx) (g> pp{t, dy). 

For (j) G Co([0,r] X M), we denote 


/ 


An{t):= dxK'^{x-y)p'^{t,dx,dy)4>{t,x)- d^K{x - y)p{t,dx,dy)4){t,x), 


/ 

Jr- 


Step 1: Convergence almost everywhere in time of An{t). 


Since dxK^{0) = 0, we have 

An{t) = [ dxK^{x - y)p^{t,dx,dy)4>{t,x) - [ dxK{x - y)p{t,dx,dy)<p{t,x), 
= ^n{t) + 

where ln{t) and II„(t) are defined by 

I„(t) := J^ (dxK^{x -y)- dxK{x - y)'j p^{t,dx,dy)(p{t,x), 

lln{t) := / dxK{x - y) {p^{t,dx,dy) - p{t,dx,dy)) (j){t,x). 

From the definition of dxK'^ in Lemma l3.11 it follows that 

^n{t) = [ {dxK^{x - y) - dxK{x - y)) p!^{t,dx,dy)(t){t,x). 

J En 

The estimate on in Lemma l3.II and (H3)| imply that 


\ln{t)\< 2 \mL^\\dxK\\L^p^{t,En), 
with /i"' and En defined in (|3.2I) . 

Let e > 0. Since the set En converges to the empty set, there exists N £N such that Vn > N, 
(3.3) p{t,En)<e. 


For all n> N, we observe that En C Ej\f, we have 

(3.4) p"-{t,En) < p^{t,EN) < - p){t,E]s[) + p{t,EN). 

From the weak convergence of p”, a = 1,2, we note that the sequence pA{t, ) converges weakly 
to p{t, •). Since the total variation of p^{t, •) is constant in n, the tight convergence is achieved. 
Then, there exists N' such that yn > N' > N 


- p\ {t,E]y) < e. 

From (j3.4p and (j3.3p . we conclude that Vn > iV' > N, 

p^{t,En)<2e. 

Hence, for all n > N' , we get 

(3.5) |I„(t)| < 2U\\L^\\dxK\\L^p^{t,En) < m\L^\\dxK\\Loo e. 

We deduce that I„(t) —>• 0. 
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Next, we show that Iln(t) tends to zero. 

Iln(t) = [ {da;K{x - y) - dxK^{x - y))(p{t,x){^/'{t,dx,dy) - fj,'^{t,dx,dy)) 

+ / dxK^{x - y)<p{t,x){iJ,'^(t,dx,dy) - fj.{t,dx,dy)), 

JR2 

where R is an integer which will be fixed later. From the construction of dxK^ in Lemma EH 
we get 


11;^ = / {dxK{x - y) - dxK^{x - y)) 4 >{t,x){n"-{t,dx,dy) - y{t,dx,dy)). 

J Eh 

Therefore, one has 

< 2\\dxK\\Lo.mL^ iy^{t,ER) + fi{t,ER)) . 

Let e > 0. Using ()3.4p . by the same token as previously, there exists N such that for all n > N, 

En) < 2 e, y{t, En) < e, 

Setting R = N, we conclude that for all n> N, 

|llUt)| <6e\\dxK\\L^\\(l)\\L^. 

For Il^(t), we notice that dxK^{x — y)(l>{t,x) is a continuous function that vanishes on the 
diagonal (x, x) and we have 


[ dxK^{x - y)(l){t,x){y'^ - y){t,dx,dy) = [ dxK^{x - y)(l){t,x){y^ - y){t,dx,dy). 
JR2 Jr2 

The tight convergence of /r” to y implies that there exists N” > 0 such that for all n > N” 


[ dxK^{x - y)4>{t,x){yJ^ - y){t,dx,dy) 
JR2 


< e. 


Therefore for all n > max{A^', iV''}, one has 

(3.6) \lln{t)\<e{l + 6 \\dxK\\L^U\\L^). 

This implies that Iln(t) converges to 0. 

Combining ()3.5p and ()3.6p . we deduce that for almost every t € [0,r], An converges to 0. 
Step 2: Lebesgue’s dominated convergence theorem 


For all t G [0,T], we have that 

|^n(t)| < ‘^\\4>\\L'^\\9xE\\L°°MaMi3. 

Since An converges almost everywhere to 0, An{t)dt converges to zero from Lebesgue’s dom¬ 
inated convergence theorem. □ 


3.2. OSL condition on the macrosocopic velocity. 


Proposition 3.3. Let T, M he positive constants and p G Sm be a positive measure such that 
Vt G [0, T], IpI (t,M) = M. Let K he such that assumption (H4) hold. Let d{p) and aL{p) he 
defined in dMl) and (EH) respectively. Then, there exists k G L^([0,T]) such that 


dxd{t,x) < K{t), dxa^{t,x) < K{t), in the sense of distributions. 
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Proof. For x, y G M, we compute: 

{a{p){t,x) - a{p){t,y)){x - y) = f {dxK{x - z) - dr,K{y - z))(x - y)p(t,dz). 

Jm 

Using the A-concavity of iF, we deduce from (12.111 

(d(p)(t,x) - d(p)(t,y))(x -y) < \{x - yf f p{t,dz) < XM{x-yf. 

Since is also A concave from the proof of Lemma 13.11 we get the one-sided Lipschitz estimate 
on a"' by the same token as for a. □ 

4. Existence and uniqueness of duality solutions 

4.1. Proof of the existence of duality solutions in Theorem 12.51 The proof is divided 
into several steps. First, we construct an approximate problem for which the existence of duality 
solutions holds. Then, we pass to the limit in the approximate problem to get the existence of 
duality solutions thanks to the weak stability of duality solutions stated in Theorem 12.31 and 
recover Equation ()2.5I1 from Lemma 13.21 Finally, we recover the bound on the second order 
moment. 


Step 1 : Existence of duality solutions for the approximate problem 


The macroscopic velocity a is replaced by an approximation a” defined in (|3.ip and the following 
system is considered: 


(4.1) dtPa + Xadx {a^{0iPi + 02P2)Pa) = 0, for 0 = 1,2. 


Since dxK"' is not Lipschitz continuous, we first consider dxK'^’^ an approximation of dxK"^ 
obtained by a convolution with a molifier. The solution pa’^ to the following equation is inves¬ 
tigated. 

(4.2) dtpl^^ + Xadx + e2p'^,nPin = 0, for a = 1, 2, 


where a”’’”* is given by 

a---(p) := f dxK^’^ix-y)p{t,dy). 

Jr 

Applying Theorem 1.1 in m gives the existence of solutions pa"* in L°°([0,T], A4b(M)) and 
|pa™'| (t,K) = |p^*| (M) = 1. Since the velocity field a"’’™' is Lipschitz, pa’^ is a duality solution. 
In addition, for f G C')T(M) we have for a = 1, 2 the following estimate: 


fL 

dt 



(t, dx)(f{x] 


dxK'^'^^ix 


y)( 6 'ip"’™(t, dy) + 02P2’”'{t, dy))p^'^{t, dx)dx(t){x). 


Then, 

(4.3) 


d 

dt 



{t, dx)(j){x) 


< lISjjArllioo (01-|-02)- 


Using (j4.3p and the density of C“(M) in C'o(M), we deduce that pa™ G Sm- Moreover, the 
equicontinuity of pa’^ in <Sxt follows from (14.3p and the density of C“(M) in Co(M). Since 
I Pa’™' I (t,K) = |p™*| (M) = 1, Ascoli Theorem gives the existence of a subsequence in m of pa™ 
which converges to a limit named p” in Sm- We pass to the limit when m tends to infinity in 
Equation (j4.2p and obtain that p” satisfies (14.ip . 


Step 2 : Extraction of a convergent subsequence of p^ and existence of duality solutions. 
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As above, there exists a subsequence of p'^ in Sm such that 

Pa Pa weakly in Sm-, for a = 1,2. 

Let us find the equation satisfied by pa in the distributional sense. Let cf) be in (^^([Ojr] x M). 
Since p^ satisfies (I4.ip in the distributional sense, we have 

/ / dt(t){t,x)p^{t,dx)dt + Xa [ [ a^iOiPi+92P2)Pait,dx)(j){t,x)dt = 0 . 

Jo Jm. Jo Jm. 

Using Lemma 13.21 we can pass to the limit in the latter equation and obtain, 

/ / dt(j){t,x)pa{t,dx)dt + Xa / / diOipi + e2P2)Pait,dx)4>{t,x)dt = 0. 

Jo Jm. Jo '/M. 

Thus Pa satishes (12.5p in the sense of distributions. From Proposition 13.31 the macroscopic 
velocity aP{p^) satisfies an uniform OSL condition. Then, by weak stability of duality solutions 
in (see Theorem 12.31 f3U. we deduce that 

dtPa + Xadx{d{p)pa) = 0, for a = 1,2, in the sense of duality in ]0,T[xR. 

Step 3 : Finite second order moment. 

From Equation p2.5p . we deduce that the first and second moments satisfy in the sense of 
distributions 

d 


a = 1 , 2 . 


\x\pa{t,dx)^ = - j sgn{x)d{p)pa{t,dx), 
\x\'^ Pa{t,dx)^ ^ 


Since € ^ 2 ( 1 ^) and d{p) is bounded from ((H2)), we deduce that the first two moments of 
Pa are finite, then Pait) € ^ 2 ( 1 ^) for t > 0. □ 

Remark 4.1. If we define the weighted center of mass of the system Xc as follows: 

Xc{t) := — [ xpi{t,dx)-\—- [ xp 2 {t,dx). 

Xi Jr X2 Jr 

We remark from straightforward computation that ^Xc = 0. Then the weighted center of mass 
is conserved for this system. 

4.2. Proof of the uniqueness of duality solutions in Theorem 12.51 Uniqueness relies 
on a stability estimate in Wasserstein distance, which is the metric endowed in P 2 (I^)- This 
Wasserstein distance dw is defined by (see e.g. [261 [27]) 

dw{v,p)= inf 1/ \y - x\^ x{dx,dy)\ 

7Gr(^,At) J 

where F(/i, v) is the set of measures on x R^ with marginals p and v, i.e., 
r(z^,/r) = {7 G 7’2(R X R),VC G C'o(R X R), f ^{yo)x{dyo,dyi) = [ ^{yo)p{dyo), 

I JR^ JR 

/ i{vi)l{dyo-,dyi) = / i{yi)iy{dyi) 

2r2 Jr 

The Wasserstein distance dw takes a more pratical form in the one-dimensional setting. Indeed, 
in one space dimension, we have (see e.g 
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where ^ and ^ are the generalised inverse of cumulative distributions of u and /i, defined 
by 

F~^{z) = inf|x G M, z^((—oo, x)) > 2 ;|, F~^{z) = inf|x € M, //((—oo, x)) > zj. 

This Wasserstein distance can be extended to the product space 7^2(1^) x 7^2 (l^)- In the case at 
hand, we define W 2 (i^,fJ.) by 

(4.4) W2(i^,f^f = 1^' \F-Hz) - F-Hz)f dz + \F-\z) - F-\z)f dz, 


where u = ^ 7^2(1^) x 7^2(1^) and F^^^, F^^ are the generalised inverse of 

cumulative distributions of Va and Ha for a = 1, 2, respectively. Using VU 2 we prove a contraction 
inequality between duality solutions of dni). 


Proposition 4.2. Let = 


, ,ini 

Ml 


and = 


, ,ini 


he in 7^2 (K)^- lUe define ji 


and u = 


duality solutions of dni) with respectively the initial data pd 


Then W 2 {p-,y) defined in (14.41) is bounded and satisfies the estimate: 



W 2 {p, v) < W 2 {r\ exp (2A(xi + X2)(0i + 02)7). 

Proof. Since (//, v) G 7^2(1^)^ x 7^2(l^)^, is bounded. For the sake of clarity in the proof, 

we denote 

Fa"--=Kl\ G-1:=F-;, for a = 1,2. 

We also omit the argument t in notations F~^{t,x) and G~^{t,x). Computing the derivative of 
W 2 {fJ^,i')‘^ with respect to time, 

dtW 2 {pi,uf = 2 [' {Ffi\x) - GfHx)) {dtFfi^x) - dtGf\x))dx 
Jo 

+ [ {Ffi^ix) - Gf^{x)) {dtFfi^{x) - dtGf^{x))dx. 

X2V1 Jo 

Straightforward and standard computations give that 

dtF-^{x) = Xad{t,F-^{x)), dtG~^ = Xad{t,G~^{x)), for a = 1,2. 

From the dehnition of d in (|2.6I) . we get 

dtFfi'^{x) = xiOi [ dxK{Ffi^{x) - z)fii{t,dz)+ X 1 O 2 [ dxK{Ffi^{x) - z)fj. 2 {t,dz). 

Jo Jo 

Setting z = Ffi^{y) in the first integral and z = Ffi^{y) in the second one yields 

dtFfi\x) = X101 [' ^iFfi\x) - Ffi\y))dy + xi02 [' ^(Ffi^x) - Ffi\y))dy. 

Jo Jo 

Similarly, we get 

dtGf\x) = xiOi ['d^{Gf\x)-Gf\y))dy + xi02 [' d^{Gf\x) - Gf\y))dy. 

Jo Jo 
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Using the oddness of dxK, we can symmetrise the terms in the right-hand side of dtF^ dtG^ 
One gets 

[' {F^Hx) - G^Hx)) {dtF^\x) - dtG^\x))dx = 

Jo 

(^d^{F,-\x) - F^\y)) - d^{Gi\x) - Gi\y))) x 
(Ff ^(x) - G^^{x) - (Ff ^(y) - G^^{y))) dxdy 

(£k(F,-‘(i) - F,-■(!,)) - OKOr'W - Oj‘(9))) (Ff'W - 0^\x))dydx. 




1 fI 


+ Xl^2 



0 JO 


Similar computations can be carried out for (F 2 ^ (t, x)—G 2 ^ ®)) (^*-^2 ^ -dtG2\t,x)). 

Finally, dtW 2 {i',y)'^ reads 

dtW 2 {u,yf = xiOi (fk{F{Hx) - F{\y)) - &(Gr'(x) - Gr'(y))) x 

(Ff ^(x) - Gj-^x) - (Ff^(y) - G5-^(y))) dxdy 
^ [WF^Hx) - F^Hy)) - WGI^Hx) - G 2 \y))) x 

{F 2 ^{x) - G 2 ^{x) - {F 2 ^{y) - G 2 ^{y))) dxdy 
+ 2x102^ ^ [d^{F-\x)-F-\y))-d:kiG^\x)-G2 \y)))x 
(Ff i(x) - F 2 "^(y) - (G5 -^(x) - G 2 ^(y))) dxdy. 

Applying inequality (j2.1l) to (j4.5p and using Young’s inequality yields 

dtW2{i^,yf < 4 xiAx ^(01 + 62) j (Ff ^(x) - G^^{x)fdx + {62 + j “ G2^{x)fdx 

By definition of 14^2 (|4.4p , we conclude 

dtW 2 iu, yf < 4A(xi + X2)(01 + e 2 )W 2 {v, yf. 

Then the result follows from Gronwall’s Lemma. □ 

Proof of uniqueness. From Proposition 14.21 it is clear that if /r*"'* = F”®, then fi = v. We 
deduce uniqueness of duality solution in Theorem 12.51 


4.3. Equivalence with gradient flow. We recall that y € AC^([0,F],^ 2 ( 1 ^) x ^ 2 ( 1 ^)) if fJ 
is locally Holder continuous of exponent 1/2 with respect to the Wasserstein distance 41^2 in 
^ 2 ( 1 ^) X F2(1K)- 


/ ^ini 

Proposition 4.3. Let assumptions of Theorem \2.5\ hold. Given p*®®® = ' ^ 


jm ; € ^2(1^) X F2(1K)- 

F2 


Let p = 


and p = 


be respectively the duality and gradient flow solution. Then, we 


have p € AG"^([0, T], ^ 2 ( 1 ^) x F 2 (K)) cLnd p = p. 
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Proof. We have that 
is bounded and 



€ L^([0, T], L^(/9i (8) /92,IR^)). This comes from the fact that dxK 



y){9ipi{t,dy) + 62 P 2 {t,dy)) 


< (^1 + 02 )- 


From Theorem 8.3.1 in m, we deduce that p G AC‘^{[Q,T],'P 2 i^) x V 2 (J^)). Since p satisfies 
()2.5p in the distributional sense, we deduce by uniqueness of such solution that ^ is a gradient 
flow solution. 

Conversely, we suppose that p is a gradient flow solution, we have that p G (^([O, T],'P 2 (K) x 
P2(1K)) and p verifies (j2.5h ~ (j2.6h . By uniqueness of the solution in Theorem 12.51 we deduce that 

p = p. □ 


5. Convergence for the kinetic model 


The convergence of the kinetic model (II.2p towards the aggregation model is analysed in this 
section. 

Proof of Theorem \ 2. 61 From the assumption + 62 ) < 1 for a = 1, 2, we obtain that To,[5] 

defined in (II.3h is positive. Since ra[5] is a bounded and Lipschitz continuous function, we get 
the global in time existence of solutions to (jl.2p and we have that /„ G C{[0,T], L°° n Lf{M.)) 
and f ff^dx < 00 . 

To prove the convergence result stated in Theorem 12.61 we consider the zeroth and first order 
moments of the distribution f^{x,v,t) introduced previously. 

Pa :=/«(1)+ /«(-!), 4 := (/a(l)-/a(-l)), for a = 1,2. 

From (jl.2p . these moments satisfy the following equations 

dtPa + dxJf, = 0, 


(5.1) 


dtJa + dxPa = -{XadxS^'Pa “ Ja), for O = 1, 2. 


From the first equation of (j5.1h . we deduce that Mt G [0,T], |p^| (t,M) = \p''a''\ (K). Therefore, 
for all t G [0,r] the sequence {Pa{t,-))e is relatively compact in — f7(W(b(]R), Cq (M)). 

Since is uniformly bounded in C'°([0,T], L^(M)), using the same token as in the proof of the 
existence, there exists pa G Sm such that 

Pa Pa weakly in Sm-, for a = 1,2. 

From the second equation of (|5.1I) . we have 

= XadxS^Pa — 2 i^tJa + ^xPa )) fo the distributional sense 
:= + RP 


We have that R^ converges weakly to zero in the sense of distributions. From Lemma 
obtains 

f f d{9ipl +92pl)4>{t,x)p%{t,dx)dt ^ f [ d{9ipi + 92P2)(t>{t,x)pa{t,dx)dt. 

Jo Jo J'K 

^ Xa^i^iPi + 92 P 2 )Pa in the sense of distributions. 


one 


We conclude that 


Passing to the limit in the first equation of (j5.1h . we deduce that pa satisfies (j2.5p in the sense 
of distributions. We use uniqueness of duality solutions to conclude the proof. □ 
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6. Numerical simulations 

This section is devoted to the numerical simulation of system (12.51) . We provide a numerical 
scheme which preserves basic properties of the system such as positivity, conservation of mass 
for each species and conservation of the weighted center of mass. Moreover, we prove the 
convergence of the numerical approximation towards the duality solution defined in Theorem l2.5l 


6.1. Numerical scheme and properties. Let us consider a cartesian grid of time step At and 
space step Ax. We denote Xj = jAx,j € Z,t"' = nAt,n G N. An approximation of paif^jXj) 
denoted is computed by using a finite volume approach where the flux Faj-i /2 given 
by the flux vector splitting method (see [E]). Assuming that {p^j) are known at time we 
compute by the scheme: 


( 6 . 1 ) 


ra,j ra,j I a,J + l/2 a,j — ll2 n £ in j ^ ^ 

—- - H- - —- -—— = 0 for a = 1, 2 and j € Z, 

At Ax 

< Kj-i /2 = 

= X] ~ {^lPl,i + ^2p2,i) , 

' i¥=j 


where (•)"'■ := max{(-),0}, (•) := min{(-),0} are respectively the positive and negative part of 

(•). Then we reconstruct 


Nt-l 

PaAxit,x)= 

n=0 jez 

where 5xj is the Dirac delta function at xj = jAx. We first verify that this scheme allows the 
conservation of the mass and of the weighted center of mass. 


Proposition 6.1. Let us consider (p)^”*, P 2 "’*) € 7^2(1^)^ such that for a = 1, 2, p)(^* = Yljsz Paj^xj ■ 
We assume that for n € N*, {Pa,j)j,n o>fe given by the numerical scheme (16.Ih . Then the con¬ 
servation of the mass of each species and of the weighted center of mass hold: 

( 6 . 2 ) VneN, YPaj" = YPaJ for a = 1,2, 


(6.3) 


Vn G N, 



iez 


h. 

X2 


Y^^iPlr 

jez 


Proof. Identity (|6.2I) can be obtained directly by summing over j G Z the first equation in (|6.1|) . 

We now show (|6.3h . Multiplying by xj the first equation in (16.Ih and summing over j G Z, 
one gets 




jGZ 


jez 




Ax 


jez 


jez 


+ for «= 1,2. 


jez jez 

Using a discrete integration by parts, one gets 


Xa 


jez 


Xa 


jez 


Xa 


jez 
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Finally, we get 

^ E ^ E pipS' = ^ E piPli + S E piPli + A* E (o^pli + ^^ph) ■ 

jez jGZ jGZ jez jGZ 

From (EH), we have that 

^ a] {e^pl^ + e^pl,) = Y, d.K{x, - X,) {9 ipIj + 02pI,) + 92Py ■ 

jGZ i^j 

Swapping indices i and j and using the oddness of dxK yields 

jez 

Then (16.,lb follows. □ 

Lemma 6.2. Let (Pi™, P 2 ™) m 7^2(1^)^ such that p^® = YljezPa j^xj with YljezPaj — ^ 

Pa,j > 0 for a = 1,2. Assuming that for n € N*, {Paj)j,n o-re given by the numerical scheme 
EH- If the following CFL condition holds 


(6.4) 


At 

cK\\L°°{dl + ^ 2 )-^ < 1 , 


Then for all n G N, p”^- > 0 and we have sup 




E \\dxK^L°°{9i + 62 )- 


Proof. This result is proved by induction. Let us assume that at time n, for all j € Z, p”^- is 


positive and sup^ 
(6.5) 


< \\dxK\\Lcx>(^9i + 02 )- From (16.11) . it follows that 




^ l“"l) Pl, + ^(p1-Pp2,-i - YS.,+.- 


Using the condition (16.4p and the fact that sup^ o” < \\dxK\\L°°{9i + 6 * 2 ), we get that 


1. Therefore P^^^ is positive as a linear combinaison of positive numbers. 

Then, recalling the expression of a” given in (16.111 . using the fact that j G Z, are positive 
and the conservation of the mass, (| 6 . 2 p yields 

<\\dxK\\LP>^{ei + e2). 


< 


0,1 


□ 


6 . 2 . Convergence of the numerical solution to the theoretical solution. In this part, 
we prove that the numerical scheme given in (|6.ip converges to the duality solution obtained in 
Theorem 12.51 


Theorem 6.3 (Convergence of the numerical scheme). Let T > 0, Ax > 0 and At > 0 such 
that (16. 4p is satisfied and denote Nt = ^- Let pj(®® G 7 ^ 2 ( 1 ^), we define 



J X . 1 


pT 


(x) dx, 


j G Z. 


Let us define pa,Ax G A4fe([0,T] x M) by 


Nt-l 

PaAx{t,x)= Y YP»d^lt’"’W+^it)^Xj{x), 

n=0 jGZ 


( 6 . 6 ) 
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where {Paj)j,n computed by ()6.1I) . 

Then, we have 

PaA^ Pet weakly in M.b{[0,T] x M) as Ax —>■ 0, 
where pa is the duality unique solution of Theorem \2. 51 with initial data p^ff^. 

Proof of Theorem \6.3i . For the initial data, it is clear that when Ax —>■ 0, we have Pa,Axil- — 
0) ^ weakly. From Lemma Wf2\. we get that for all j € Z,re G N, values of are positive. 

Step 1 : Extraction of a convergent subsequence 
Equation (16.21) implies that the total variation of pa,Ax is fixed and independant of Ax. 

|Pa.Ax|([0,r]xR)=T^pi™. 

jez 

Therefore, there exists a subsequence of paAx that converges weakly to pa G Adft([0,T] x M). 


Step 2 : Modihed equation satisfied by Pa,Ax 
Let be (/) G C“((0,T) x M). From the dehnition of pa,Ax hi (| 6 . 6 I) . we have 


n Nt -1 

dtPa,Axi4^ P — / Pa,Axl^t4^ — ^ ^ ^ ^ Pa,i I 

j[o,r]xR Q -py Jttt 


dt(pixj,t)dt, 

n=0 j'eZ 

Here and below we use < •, ■ > to denote the dual product in the sense of distributions. Discrete 
integration by parts yields 

Nt-l Nt 

n=0 jez n=ljez 

where we use the notation := (j){xj,t"'). Using (|6.5I) and applying transformations to indices 
yields 

At At 

< a,p^Ar.<P>= x; E T-ippPlPlti -+ w E E(»”)^<.i('^r‘ - 


n=0 jgZ n=0 jgZ 

Taylor expansions gives the existence of in (xj,Xj_|_i) and in (xj_i,Xj) such that 




l^n+l _ j^n+1 


in+1 


Putting together, one obtains 


, - Axd^epT^ + 

3 ^^3 2 


Nt-l 


< dtpaAxA>= At Y^lPaJ^xl^Y + Riit^X,At), 

n=0 jez 

where R]^{Ax, At) is given by 

Nt-l / t A _„',2 


At 


RA^At) :=^ E 

n=0 jGZ 

At 

E T,(pjrpz 

n=0 jez 


(Ax)" 


■dxxfiiC^.'t 


3 +n+l' 
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From (j6.6p and the definition of a in (|2.6p . we have 

Nt-l 

+ ^ 2 / 52 ,Ax) = E r>.+ -L[it)6xjix). 


n=0 j£Z 


where are defined in (jd.ip . We get that 


j-n +1 


< Ax + 6*2/02,Aa:)da,Ax, >= “ V] Pa,j / dx(f>{Xj,t)dt. 

n=0 jeZ 

From the Taylor expansion of dx4>{xj,t): 

dx4>{xj,t) = dx(l){xj,t'^^^) + {t- 

with G one sees that 

Nt-l 

< a(6*lPl,Ax + 02P2,Ax)Pa,Ax,dx(i) >= “At ^ ^ Oj p'^jdxCp]^^ + RI{Ax, At), 

n=0 jgZ 

where R^{Ax,At) is defined as follows: 

Nt-l 

Rl{Ax, At) = rf )(it. 

n=o jez 

The modified equation satisfied by pa,Ax in the distributional sense writes: 

/ / / 0 a,Ax 5 t(/)(f, x) + / / 0(01^1,Ax + 6*2/02,Ax)/ 0 a,Ax 5 x'/> = -Ra( Ax, At) + (Ax, At). 

</0 i/M. </0 i/M 

From Lemma 16.21 we deduce that the terms R)^ and R^ satisfy the estimates: 


\Ri\ <CTAx\\dxxcji\\L^, 


\Rl\ < CTAxWdtxHL^ 


where C stands for a nonnegative constant. Passing to the limit and using the technical Lemma 
Mi we conclude that the limit pa satisfies (|2.5I) in the distributional sense with the expression 
p2.6p for the velocity. By uniqueness result in Theorem 12.51 we deduce that pa is the unique 
duality solution of (II.ip . □ 

6.3. Dynamics of aggregates and numerical simulations. In this part, we carry out sim¬ 
ulations of Equation (|2.5p obtained thanks to scheme dEU). Before numerically simulating the 
hydrodynamic behavior of the chemotaxis model, we first clarify the aggregate dynamics of this 
model, especially on the synchronising dynamics between aggregates of different species. 

For the sake of simplicity, we choose 6 i = 62 = ^ and K = in (j2.5h . which corresponds to 

the particular case of bacterial chemotaxis (see (11.71) 1. To illustrate the synchronising dynamics 
of the aggregates for (12.5p . we consider the initial data given by sums of aggregates 

k k 

and look for a solution in the form 

^ ^ (t) 1 ^ ^ (i) ‘ 

k k 

We denote by ui and U 2 antiderivatives of pi and p 2 , respectively. Then the equation p2.5l) reads 
(6.7) dtUa +Xadpa = 0, a = 1,2, 
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in the sense of distributions. Direct computation shows that 

api = y^^pkjpedxKjxk - xi) + i^AKixk - ye))Sxk, 

k,£ 

ap2 = ^ Uk{pedxK{yk - xi) + - yi))Sy^- 

k,l 


Injecting these expressions into equation (16.71) . the positions x^ and y^ satisfy the system of 
DDEs 


x'kit) = Xi '^iP£dxK{xk - Xi) + UidxK{xk - yi)), 
£ 


y'kit) = X2 '^{p£dxK{yk - Xi) + UidxK{yk - yi)). 
£ 


We recover the same system for particle solutions as in DiFrancesco and Fagioli [n] for two 
species. See also similar aggregate dynamics for single species in um- In the case of one 
single species, the system of ODEs is determinant before any collision of aggregates, and after 
each collision, one can always ‘restart’ the particle system till final collapse of all aggregates. 
However, the case of collisions between particles of different species is more complex, since it 
does not necessarily imply whether the particles of different species will synchronise or not 
after colliding. In fact, as observed in the following simulations, both ‘synchronising’(colliding 
particles of different species staying together) and ‘non-synchronising’ cases can occur, and the 
transitions between the synchronising types may happen, depending on the weighted attraction 
of other aggregates acting on them. 


For illustration, we assume that two points of different species collide at a time to- For 
instance, take Xk{to) = yki^o) for some k, then at this time to we have 


(6.8) x'kito) = xilkito), y'kito) = X2'yk{to), Ik = '^iP£dxK{xk-Xi)+UidxK{xk-yi)). 

£,k 


Note that Xk characterises external weighted attraction on Uk and pk: depending on chemo- 
sensitivities, distances to other aggregates and the masses of all other aggregates. 

Thus if xi 7 ^ X 2 the velocity of species 1 and 2 is not the same at this time to- However, 
with the special case at hand, K{x) = we have dxK{xk — yk) —>■ \ when Xk — yk] and 

dxK{xk — yk) when Xk ^ yk- We deduce that when Xk < yk and Xk yk we have 

{yk - XkYit) = -^{xiJ^k + X2Pk) + {X2 - Xi)lk{t). 

Obviously, in this case, particles pk and Vk stay together if {yk — Xk)'{t) < 0 . On the other hand, 
when yk < Xk and yk —?> Xk we have 

{xk - 2/fc)'(t) = -^(xWfc + X2Pk) + (Xi - X2)lk{t)- 

In this case, particles pk and Uk stay together when {xk — yk)'{t) < 0 . Finally, to keep Xk{t) = 
yk{t) for t > to, we need the condition 

(6-9) Kxi - X2)lk{t)\ < ^(xiz^fc + X2Pk), 

where 7 fc(t) is defined in (16.81) . This relation characterises the weighted attraction of other 
aggregates acting on them. If the external weighted attraction on Uk and pk (the left hand side 
of (|6.9p ) is small, they will stay together. When the external weighted attraction is big, the 
attraction between I'k and pk is relatively weak and they will move separately, the one with 
bigger motility will move faster than the other. 


SYNCHRONISING AND NON-SYNCHRONISING DYNAMICS 


19 


We call (j6.9p the synchronising condition for and ■ Similarly, we can get the synchronising 
condition for any //j and Uj, If more than two aggregates collide simultaneously, we can 

simply replace them by two aggregates of each species, each aggregate accumulating the total 
mass of each species. 

In conclusion, according to the dynamics defined above, we can see that the initial aggregates 
will collapse such that they eventually form a single aggregate of the two species. The final 
aggregate can not separate, which is similar but illustrate more complex behaviour as one 
species case discussed in [18]. Now we give some numerical examples showing “synchronising”, 
“non-synchronising”, transitions between “synchronising” and “non-synchronising” dynamical 
behaviours for the hydrodynamic model ( 12 .5p . 

Example 1: Synchronising dynamics. Take the chemosensitivity constants xi = 10) 
X 2 = 1 in ([ 22 ]), and consider initial data 

^0 _ ^g-5000(x-|-0.5)2 _|_ 2g-5000(a;-0.5)2 ^0 _ 2g-5000(a:-l-0.15)2 

It corresponds to small bumps located at position xi(0) = —0.5, yi(0) = —0.15, X 2 ( 0 ) = 0.5, 
with mass /ii = 4mo, //2 = 2mo, fi = 2mo, where mo = dx. Figure [ 6 T] displays 



t= 1.089611e+00 



t= 1.231734e+00 



t= 1.373857e+00 

6 I-.- 1 —.-.-1 

5 ■ 

4 • 

3 ■ 

2 • 

1 • 

0 I - . -L—.- . - 

-1 -0.5 0 0.5 1 


Figure 6.1. Example 1. Snapshots of pi (red solid line) and p 2 (blue dashdot). 

The evolution shows the synchronising dynamics after first collision. 

numerical results obtained thanks to the scheme dEI]) defined above. We first observe the fast 
blow-up with the formation of Dirac deltas. Then, the numerical simulation shows that pi and 
vi collapse for the first time at ti « 0.947, with xi(ti) = yi(ti) ~ —0.18, and X 2 (ti) ~ 0.12. We 
check the “synchronising condition” (|6.9p : 

LHS = Kxi - X2)7i(ii)l = (10 - 1) X 2 X = 9 e-(" 2 -i) < 9 , Vxi,X 2 , 

RHS = i(xwi + X2f^i) = ^(10 X 2 + 1 X 4) = 12. 

Thus the “synchronising condition” (16.9p is always satished, then they will move together af¬ 
terwards till final collapse with p 2 - This evolutionary dynamics is shown in Figure [ 6 Tj The 
numerical result confirms the synchronising dynamics of the aggregates. 
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Example 2: Non-synchronising dynamics. Take the chemosensitivity constants xi = 10; 
X 2 = 1 in (|2.5p . and consider initial data 

^0 _ 2g-5000(3;+0.5)2 ^^-5000{x-0.5f ^0 _ 2g-5000(a;+0.15)2 

It corresponds to small bumps located at position xi(0) = —0.5, yi{0) = —0.15, X2(0) = 0.5, 
with mass = 2mo, ^2 = 4mo, = 2mo, where mo = Jg dx. The numerical simulation 

in Figure [62] shows that yi and i^i collapse for the first time at ti ss 0.9, xi(ti) = yi(ti) ~ —0.15, 
and X 2 (ti) ~ 0.25. Direct computation shows that 

LHS = Kxi - X2)7i(ti)l = (10 - 1) X 4 X « 18e-°-^ 12.066, 

RHS = ^(xii^i + X 2 yi) = ^(10 X 2 + 1 X 2) = 11, 

thus the “synchronising condition” (16.9p is not satished, then they will change their order after 
intersection and travel separately. The simulation shows will collapse with //2 at time ^2 ~ 
1.61, and finally collapse with i^i at time ~ 1.85. This dynamics is shown in Figure (621 The 
numerical result confirms the non-synchronising dynamics of the aggregates. 
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Figure 6.2. Example 2. Snapshots of pi (red solid line) and p 2 (blue dashdot). 
The evolution shows the non-synchronising dynamics after hrst collision. 


Example 3: Transition from synchronising to non-synchronising dynamics. Take 
the chemosensitivity constants xi = 10, X 2 = 1 in (12.5p . and slightly modify the initial data of 
Example 2 to 

^0 _ 2g-5000(x+0.5)2 ^g-5000(a;-0.5)2 ^0 _ 2g-5000(x+0.3)2 

It corresponds to small bumps located at position xi(0) = —0.5, 2/i(0) = —0.3, X2(0) = 0.5, 
with mass = 2mo, P 2 = 4mo, ui = 2mQ. The numerical simulation displayed in Figure iQl 
shows that pi and ui collapse for the first time at ti ps 0.47 with a:i(ti) = yi(ti) ~ —0.29, and 
X 2 {ti) ~ 0.39. Direct computation shows that 

LHS = Kxi - X2)7i(ii)l = (10 - 1) X 4 X « 18e-°-®® « 9.1191, 

RHS = ^(xiui -b X2Pi) = ^(10 X 2 -b 1 X 2) = 11, 

thus the “synchronising condition” ()6.9p is satisfied, then they will move together toward p 2 - 
The interesting phenomenon is that, as their distance to p 2 is decreasing, the LHS of the 
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“synchronising condition” ()6.9p is increasing and finally greater than the RHS. The simulation 
shows that, at t 2 ~ 1.04, xi{t 2 ) = 1 / 1 ( 12 ) ~ —0.26, and X 2 {t 2 ) ~ 0.23, then 

LHS = Kxi - X2)7i(li)l = (10 - 1) X 4 X ^ ^ 

then after this time t 2 -> the interaction type has been changed: the “synchronising condition” 
f|6.9p is no longer satisfied, then they will travel separately. Further simulation shows that 
collapses with /X 2 at time ~ 2.037, and finally collapse with ui at time ~ 2.32. The full 
dynamics is shown in Figure [6131 The numerical result shows the transition from synchronising 
to non-synchronising dynamics of the aggregates. 
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Figure 6.3. Example 3. Snapshots of pi (red solid line) and p 2 (blue dashdot). 

From time to = 0 to ti ^ 0.47, pi moves toward ui. From time ti ^ 0.47 to 
t 2 ^ 1.04, Pi and ui travel together. The synchronising type changed at ^ 2 - 
After time ^ 2 , pi overtakes ui and collapse with p 2 at ts ~ 2.037, and finally all 
the aggregates collapse at ^4 2.32. The evolution shows the transition from 

synchronising to non-synchronising dynamics. 

Example 4: More complex transition. Take the chemosensitivity constants xi — 10? 

X2 = 1 in (ESD, and consider initial data 

pO ^ 3g-5000(x-h0.8)2 ^ 5g-5000(a;-r0.02)2 = 3 5g-5000(x-0.02)2 g ^^-5000{x-0.5)^ 

It corresponds to small bumps located at position xi(0) = —0.8, * 2 ( 0 ) = —0.02, yi(0) = 0.02, 

2 / 2 ( 0 ) = 0.5, with mass pi = 3mo, fJ -2 = 1.5mo, I'l = 3.5mo, 1^2 = 8.5mo. The snapshots 
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of Pi and p 2 are shown in Figure 16.41 We observe that p ,2 and vi meet for the first time at 
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Figure 6.4. Example 4. Snapshots of pi (red solid line) and p 2 (blue dashdot). 


ti ~ 0.0459 and satisfy non-synchronising condition so they separate after ti. See the snapshot 
at t 2 ~ 0.4288 for evidence. They meet for the second time at ~ 0.9494 but the synchronising 
type has been changed: now they satisfy synchronising condition thus they travel together 
afterwards. At time t 4 1.04, pi catches p 2 and ui. Now we treat them as {pi -|- P 2 ) and ui: 
they satisfy non-synchronising condition and separate, see snapshot at ts 1.256 for evidence. 
At time « 1.684, {pi + P 2 ) collapse with 1 ^ 2 , satisfying synchronising condition and staying 
together till final collapse with ui at tj ~ 2.756. The illustration shows the complex changing 
of interaction types for the aggregate dynamics of two species chemotaxis model. 
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